# TABLE B.4
rm(list=ls())

(WD <- getwd())
if (!is.null(WD)) setwd(WD)

indata0      = paste0(WD,"/_individual_data/_spfrawdata/","", collapse = NULL)
indata1      = paste0(WD,"/_average_data/_spfrawdata/","", collapse = NULL)

library(haven)
library(AER)
library(sandwich)
library(lmtest)
library(pracma)
library(stargazer)
library(plm)
library(pracma)
library(DataCombine)
library(jtools)
library(plyr)
library(ggplot2)
library(tidyverse)
library(DescTools)

beta_coef_save = rep(NA, 4)
std_coef_save  = rep(NA,4)
nn_save        = rep(NA,4)
r2_save        = rep(NA,4)
f_save         = rep(NA,4)

##[1] Baseline and Winsorized 
load(paste(indata0,"us_spf_pgdp_ind.Rda",sep=""))
data$ID         = data$id
data            = data[data$qdate>=1970.00,]
data            = data[data$qdate<=2020.00,]
data            = ddply(data,.(qdate),transform, cons = mean(f2,na.rm = TRUE))

delta            = plm(fc_err ~ cons, data=data,index=c("ID", "qdate"), model="within")
beta             = coefficients(delta)
obs              = nobs(delta)
stderror         = sqrt(diag(vcovDC(delta, type = "HC1")))
r2               = r.squared(delta)
ff               = summary(delta)$fstat

beta_coef_save[1] = as.numeric(beta[1]) 
std_coef_save[1]  = as.numeric(stderror[1]) 
nn_save[1]        = as.numeric(obs[1])
r2_save[1]        = as.numeric(r2[1])
f_save[1]         = as.numeric(ff[2])

data$fc_err_win = Winsorize(data$fc_err, probs = c(0.01, 0.99), na.rm = TRUE)

delta_win        = plm(fc_err_win ~ cons, data=data,index=c("ID", "qdate"), model="within")
beta             = coefficients(delta_win)
obs              = nobs(delta_win)
stderror         = sqrt(diag(vcovDC(delta_win, type = "HC1")))
r2               = r.squared(delta_win)
ff               = summary(delta_win)$fstat

beta_coef_save[4] = as.numeric(beta[1]) 
std_coef_save[4]  = as.numeric(stderror[1]) 
nn_save[4]        = as.numeric(obs[1])
r2_save[4]        = as.numeric(r2[1])
f_save[4]         = as.numeric(ff[2])

##[2] Consensus -i 
data            = ddply(data,.(qdate),transform,cons_all = mean(f2,na.rm = TRUE))

id_all          = unique(data$ID)
iter            = 2
data_old        = data
for (ii in id_all){
  print(ii)
  data_tmp    = data_old[!data_old$ID==ii,]
  data_tmp    = aggregate(. ~ qdate, data_tmp, mean)
  data_tmp    = subset(data_tmp, select = c(qdate, f2))
  names(data_tmp)[iter] = paste("cons",ii, sep = "")
  
  data        = merge(data,data_tmp, by = c("qdate"))
}

beta_ind_save   =  rep(NA, length(id_all))
nobs_win        =  10
iter            =  0
for (ii in id_all){
  print(ii)
  iter         = iter + 1
  data_tmp     = data[data$ID==ii,]
  nobs         = nrow(data_tmp)
  tmp1         = paste("cons",ii, sep = "")
  data_tmp$tmp = data_tmp[[tmp1]]
  if (nobs>= nobs_win) {
    tryCatch({plm_tmp    = lm(fc_err ~ tmp, data=data_tmp)
    beta                 = coefficients(plm_tmp)
    beta_ind_save[iter]  = as.numeric(beta[2])}, error = function(e){})
  }
}
median_delta       = median(beta_ind_save, na.rm=TRUE)
beta_coef_save[2]  = median_delta 


##[3] Previous Consensus Deviation
data$cons_dev      = data$f1-data$cons

delta_clemens       = plm(fc_err ~ cons_dev, data=data,index=c("ID", "qdate"), model="within")
beta                = coefficients(delta_clemens)
obs                 = nobs(delta_clemens)
stderror            = sqrt(diag(vcovDC(delta_clemens, type = "HC1")))
r2                  = r.squared(delta_clemens)
ff                  = summary(delta_clemens)$fstat

beta_coef_save[3] = as.numeric(beta[1]) 
std_coef_save[3]  = as.numeric(stderror[1]) 
nn_save[3]        = as.numeric(obs[1])
r2_save[3]        = as.numeric(r2[1])
f_save[3]         = as.numeric(ff[2])

##[4] Table
all_tables            =  data.frame(beta_coef_save, std_coef_save, nn_save,f_save,r2_save)
rownames(all_tables)  = c("Prev. Consensus", "Prev. Consensus -i","Prev. Consensus Dev.", "Prev. Consensus Wins.")
colnames(all_tables) = c("Coef.", "Std.error", "#Obs.", "F-stat","R2")

write.table(format(all_tables,digits=3), file = "_tables/table_b4.txt", sep = "\t"
            
            )